Mapping molecular models to continuum theories for partially 

miscible fluids 



Colin Denniston 1,2 and Mark 0. Robbins 1 
1 Department of Physics and Astronomy, 
The Johns Hopkins University, Baltimore, Maryland 21218, USA and 
2 Department of Applied Mathematics, 
The University of Western Ontario, 
London, Ontario N6A 5B8, Canada 
(Dated: February 2, 2008) 

Abstract 

We map molecular dynamics simulations of fluid-fluid interfaces onto mesoscale continuum the- 
ories for partially miscible fluids. Unlike most previous work, we examine not only the interface 
order parameter and density profiles, but also the stress. This allows a complete mapping from 
the length scales of molecular dynamics simulations onto a mesoscale model suitable for a lattice 
Boltzmann or other mesoscale simulation method. Typical assumptions of mesoscale models, such 
as incompressibility, are found to fail at the interface, and this has a significant impact on the 
surface tension. Spurious velocities, found in a number of discrete models of curved interfaces, are 
found to be minimized when the parameters of the mesoscopic model are made consistent with 
molecular dynamics results. An improved mesoscale model is given and demonstrated to produce 
results consistent with molecular dynamics simulations for interfaces with widths down to near 
molecular size. 

PACS numbers: 68.05.-n,47.11.+j,83.10.Gr,83.10.Mj 
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I. INTRODUCTION 

Mesoscale continuum models are increasingly popular for studies of complex fluids 
Rather than using a constitutive equation for the local stress based purely on some function 
of strain, these models incorporate a dependence on the internal microstructure by including 
the evolution of a local order parameter. They have met with success in describing bulk 
properties of some materials, such as shear thinning in liquid crystals [3|, shear banding 
flows y, and phase ordering in binary fluids [5|. However, their treatment of interfacial 
stresses has not been tested for consistency with large scale molecular simulations. This is 
an important omission since the detailed interfacial behavior can have a dramatic influence 
on macroscopic flows. 

One example where molecular scale interfacial properties are important is in pinchoff of 
fluid drops. In cases where a fluid drop breaks up due to some external force, there can be a 
cascade of instabilities down to microscopic length scales jsj. How such instabilities are cutoff 
is a question of active research |2[ with practical applications to coatings of micro-particles 

II . Another example is dynamic wetting jjj. When a liquid-liquid interface intersects 
a stationary solid boundary it makes a well defined angle with the solid known as the 
contact angle. When the solid is moving, the dynamic contact angle 9d is a function of the 
wall velocity. There is significant interest in reproducing this velocity dependence within 
mesoscopic models [10]. However, 6 d is affected by details of the fluid-fluid interface and 
the solid-fluid interface that are currently unknown. Quantitatively reproducing these effects 
requires a detailed examination of the microscopic structure of the interfaces near the contact 



point and how this information can be mapped to the scale of continuum models [11 1. 

The difficulty in selecting an appropriate model is that many different parameters in the 
mesoscale model influence the value of some macroscopic property like the surface tension. 
As a complete set of parameters is unavailable, researchers are led to pick parameter sets 
primarily for convenience, rather than based on some underlying knowledge base. This can 
lead to significant uncertainties as to the length and time scales for which these models are 
applicable. Our aim in this paper is to resolve some of this uncertainty in parameter choice, 
at least for the simple fluid model we examine here. This should then serve as a guide for 
more complex fluids. 

In this paper we examine a binary mixture of simple fluids. The time and distance scales 
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involved in the hydrodynamic flow and diffusion of such fluids are accessible to molecular 
dynamics simulations. This allows us to map out the parameters of the mesoscopic model 
from the molecular simulations. In addition, we can test this mapping by comparing its 
predictions to simulations for a range of situations not explicitly used in the fits. As methods 
to map parameters from microscopic to macroscopic models are not well established, we feel 
it is essential to use simple models that allow extensive testing of the mapping. 

The final mesoscale model matches changes in the local stress, as well as order parameter 
and density profiles, through interfaces in the system. Non-local terms in the free energy, 
in this case gradients of order parameter and density, are essential for reproducing the 
observed microscopic stress at interfaces. The results show that many common assumptions 
are invalid. For example, many models neglect density variations because the bulk fluids are 
essentially incompressible. Despite this, we find that density variations at the interface still 
have a significant effect on the interfacial tension. Perhaps more surprising is that some of 
the elastic constants multiplying gradient terms are negative: The system is stabilized by 
atomic discreteness at short scales. 

In the next section we outline our molecular and mesoscale models, along with the meth- 
ods we use to simulate them. In Section ITTT1 we describe the molecular dynamics characteri- 
zation of both the bulk phases and interfaces of the binary fluid. Fits to various free energy 
functionals are described in Section IIV1 A comparison to some simulations of situations 
not used in the fitting procedure is given in Section We conclude with a summary and 
discussion of implications for further work. 



II. SIMULATIONS AND MODELS 
A. Molecular Scale Model 



As we wish to explore a wide range of static and dynamic properties in our model, it 
is essential that it is kept relatively simple. We will use a model similar to that used by a 
number of researchers to examine critical properties of fluid mixtures [3], capillary waves 
[3], and interfacial slip [3, llfij . The model consists of a mixture of two types of 
molecules, labeled 1 and 2. The interactions between atoms of type % and j separated by a 
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distance r are modeled using a Lennard- Jones (LJ) potential, 

Vv(r) = 4ev [{a^/rf 2 - fo/r) 8 ] , (1) 

where specifies the interaction energy and cry the interaction length. Unless specified, the 
force is truncated at the minimum of the potential r c = 2 1 ' 6 t7jj to improve computational 
efficiency. As a result, the interactions are strictly repulsive. For simplicity, all molecules 
have the same mass m. For the cases studied here we take an = cx 22 = cr 12 = a. Two 
molecules of the same type interact with the energy scale en = e 22 = e- An extra repulsion 
is added between unlike molecules ei 2 = e 2i = e(l + e*) to control miscibility. 

The degree of phase separation is quantified by the order parameter $ = p\ — p 2 , where 
Pi and p 2 are the local densities of species 1 and 2. Schematics of the phase diagram as 
a function of density p, temperature T, and miscibility parameter e* are shown in Fig. 
For sufficiently high densities the two fluids change from completely miscible to completely 
immiscible as e* increases (Fig. H^a)) or T decreases (Fig. H^b)). For sufficiently high e* 
one can observe coexistence of two partially miscible states by adjusting the temperature or 
density (Fig. [Tfc) and (d)). 

The object is to map the molecular model to a continuum mean field theory. Thus, we 
wish to avoid close proximity to a critical point where the behavior is dominated by large 
scale fluctuations. However, if one is too far from the critical point the interface may be too 
sharp to be described by mesoscopic models with a continuously varying order parameter. 
While not obvious a priori, we find that there is a significant range of parameters where 
mean field theory can be used. The constraints would be even more relaxed in a polymer 
mixture where, due to the polymer length, the interfaces are broader and the system is more 
mean-field like. However, the convective-diffusive hydrodynamics of polymer molecules of 
any significant length are currently inaccessible to molecular dynamics time scales. 

Periodic boundary conditions are applied in all three directions of the simulation cell. In 
systems with interfaces, we take the period in the x-direction L x to be 3 — 6 times larger than 
the periods in the other directions, which are usually equal. Most results are reported for 
systems with 16384 molecules that are roughly evenly divided between type-1 and type-2. 
The total mean density was varied between p = 0.83 m/a 3 and p = 0.925 m/a 3 . For the 
system with p = 0.85 m/a 3 , the system size was L y = L z = 16. lex and L x = 74.36a. A 
typical system configuration is shown in Fig. |21 
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FIG. 1: Schematic slices through the phase diagram at (a) constant temperature, (b) constant e*, 
(c) constant e* and p, and (d) constant T and e*. In (c) and (d) the horizontal lines tie coexisting 
phases and x marks a critical point. 




FIG. 2: A liquid-liquid interface between Lennard-Jones fluids at density p = 0.85m/ cr 3 . System 
size is L x = 74.36 a and L y = L z = 16.1 a, and there are periodic boundary conditions in all three 
directions. The temperature is ksT/e = 1.1 and e* = 5. 

To speed up the equilibration time, we supplement the molecular dynamics moves with 
a Monte Carlo procedure. Every 500 steps, molecules are relabeled (from 1 to 2 or visa 
versa) according to the Metropolis transition rule. This is a procedure first used for lattice 
simulations [12J and later extended to molecular dynamics simulations of polymer mixtures 
j"!^ . After equilibration, the Monte Carlo routines are turned off. 

A characteristic time scale is given by r = o{mj t) 1 ! 2 . The molecular dynamics simula- 
tions use a time step dt = 0.007r. Typically, we run the simulation for 5 x 10 5 dt with the 
Monte Carlo routines turned on. We wait a further 5 x 10 5 dt to ensure that any convective 
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flow caused by the relaxation of the initial (nonequilibrium) state has died away. Data is 
then collected and averaged over the next 10 7 dt. 

The motion in the ^/-direction, which is always perpendicular to any interfaces in the 
system, is coupled to a Langevin thermostat [13, [l8|- The drag term in the thermostat 
has a damping rate of 0.5r _1 . This corresponds to under-damped motion so that even in 
the ^/-direction the motion is dominated by inertia. The equations of motion are integrated 
using a velocity- Verlet algorithm jisl (l9^. 

The local pressure tensor in a molecular dynamics simulation has the following form ^| 

PM = M(r)-J](^^ jdr' p 8{v-v')) 

■ r ij JCa 

= (fW a v p )(v) + J2{fija [ rfr^(r-r')), (2) 

i,j>i JC ii 

where (••■) denotes a thermal average, Cij is a contour joining atoms i and j separated 
by the vector r^, and Greek indices indicate components along the x, y, and z directions. 
In what follows, we adopt the summation convention for repeated indices. An apparent 
ambiguity in Eq. (j2J) arises from the fact that any contour would appear to be acceptable, in 
that it satisfies microscopic momentum balance. However, requiring the microscopic stress 
to conform to the symmetry and transformation properties obeyed by the macroscopic stress 
makes the contour choice unique. The appropriate contour is just the straight line between 
the two atoms, a choice originally proposed by Irving and Kirkwood |2j|. This is true 
whether or not you impose the additional requirements on the scale of the averages in Eq. 



|21| or on the corresponding microscopic instantaneous many-body variable inside 



the averages 



221. 



Calculating and binning the stress to get local information can be computationally inten- 
sive as it needs to be done for every force pair in the system. We developed a new method 
to minimize this effort. The contour line defined by ry is divided into discrete segments 
(typically four). The contribution from each segment of the contour is then binned into a 
spatial grid every four time steps. We choose the number of segments to be high enough to 
ensure that in steady-state systems the component of the pressure tensor perpendicular to 
fluid/fluid interfaces is constant, as required by conservation of momentum. 
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B. Mesoscale Model 



Continuum hydrodynamics is based on conservation laws and the assumption of local 
equilibrium. The assumption of local equilibrium requires the existence of a free energy 
density. Further, this free energy density must be expressible as a functional of the density 
p, temperature T and order parameter $. Once determined, this free energy density can 
be used to calculate the stress tensor and chemical potential required in a full mesoscopic 
continuum theory. 

A number of mesoscale theories are available. Density functional theories based on a 
weighted density approximation can reproduce the qualitative shape of order parameter and 
density profiles through static interfaces, although the stress profile does not appear to have 
been studied 24 , 2^] . However this approach is computationally challenging because it 
requires the use of integro-differential equations. Our interest is in models that are easily 
generalized to nonequilibrium situations so we restrict ourselves to the simplest "non-local" 
theory, namely one involving local densities and their derivatives. Implementations of such 
theories come under many names in the literature. Two of the most common are lattice 
Bo ltzm aen —ees and d lffuS e interface n ydrod „ Q. 

To lowest order in gradients of the total density p and the order parameter $, the free 
energy functional can be written as 

T = k B T j dr!^P + ^K p (\7p) 2 + 

^V$) 2 + i^Vp-V$j, (3) 

where T is the temperature, ks is Boltzmann's constant, and the local bulk free energy ijj and 
elastic constants K may be functions of p, $ and T. If the two phases have the same density 
and elastic constants ("symmetric" case), then the free energy density must be an even 
function of $. This requires K p q> to be identically zero or an odd function of $. Mappings 
of the elastic constants used here to those of a few other common parameterizations of the 
free energy are given in Appendix 1X1 

Conservation of the number Mi of particles of each type imposes the constraints 

Mi + M 2 = J drp, (4) 

Mi -M 2 = ! dr§ (5) 



Thus at equilibrium we wish to minimize the Lagrangian 



L = T + p p k B T \Mx + M 2 - J drp^j 



+ p^k B T Mi - M 



(6) 



where p p and p<$> are Lagrange multipliers, and the (constant) factors of k B T are included 
for later convenience. The Euler-Lagrange equations give 



H* = |^ - A $ V 2 $ - A>V 2 p- 



<9$ 



(7) 
(8) 



In deriving Eq.fjHJ) we have ignored terms arising from variations of K p and Kq, with p and 
$. Only the variation of K p <$> with $ is retained as it must be an odd function of $. The 
equations obtained upon relaxing these assumptions and the justification for ignoring the 
extra terms that result are given in Appendix [XJ 
The Lagrangian density C, 



= $ + l -K p {Vp) 2 + \k^) 2 + A>V$ ■ Vp 



-p p p - 



(9) 



does not contain any explicit dependence on spatial coordinates. As a result, Noether's The- 
orem can be used to determine the pressure/stress tensor P. It implies that conservation 
of momentum takes the form 

V-P = (10) 

where 

P aP = -C8 a p + d aP —^ + d a ^^^. (11) 



As this is related to conservation of momentum, we can associate P with a nondissipative 
stress, or pressure tensor. Substituting C and the Euler-Lagrange equation for p p gives 



Pa/3 



Po 



S a/3 + K p 



(d a p)(d p p) - -(d 1 p) 2 5 a p 



k B T 
+A $ 

+K p *[(d a $)(dpp) + (d a p)(dp$) 
-(d y $)(d y p)6 af ,] , 



(12) 
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where the hydrostatic pressure po (trace of \P a p) is 

p = k B T(pp p + $ps,-ip), (13) 

and the expressions for //$ and \i p are given in Eqs. (JJJ) and (JSJ). The relations for the 
pressure tensor P a p and the chemical potential /i$ will allow a direct comparison between 
the molecular and mesoscopic scale. 



C. Lattice Boltzmann Algorithm 

Mesoscale schemes, of which lattice Boltzmann algorithms are an example, have proved 
very successful in simulations of complex fluids 0. Lattice Boltzmann algorithms can be 
viewed as a slightly unusual discretization of the equations of motion or as a lattice version 
of a simplified Boltzmann equation. In the free energy implementation, equilibrium prop- 
erties of the fluid are naturally incorporated into the algorithm using ideas from statistical 
mechanics l26||. We use a nine velocity model on a square lattice. Minor problems related 

n n 

to Galilean invariance |26J are removed via a correction term in the second moment |28| . 

The lattice Boltzmann scheme simulates the full dynamical equations of convective- 
diffusive hydrodynamics. These include the continuity equation 

dtp + d a pu a = 0, (14) 

where u is the fluid velocity and Greek indices indicate directions. Conservation of momen- 
tum takes the form 



d t pu a + dppu a u p = -d p P a p 

(15) 



l - 3^)5 a/3 <9 7 w 7 + daUp + d p u a 



where the parameter r/ sets the viscosity r\. Finally, there is the convection-diffusion equa- 
tion, 

d t <$> + d a $u a = r g (Yv 2 /i* - d(s (®d a P a pjY (16) 

The second term on the right, describing a flow- induced diffusion |29j], is a common feature 
of lattice Boltzmann schemes and is usually negligibly small. We set the algorithmic time 
constant r g to unity so that V is the diffusion constant. Viscous stresses associated with 
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velocity gradients are not relevant here as all velocities are zero. However, we will make use 
of these terms in future works. 

To compare the results of the lattice Boltzmann simulations to molecular dynamics, it is 
useful to compare the units used in the simulations. In molecular dynamics we use Lennard- 
Jones units, that is m, a, and r. Time is discrete with steps of dt = 0.007r. The size of the 
time step is limited by the requirement that molecular motion on the steep repulsive part 
of the intermolecular potential must be resolved. In the lattice Boltzmann simulations a 
discretization of both space Ax and time At is required. We will choose Ax ~ a primarily 
based on the requirement of resolving the interface width. One could, in principle, use a 
much coarser lattice in the bulk regions. Stability of our lattice Boltzmann scheme requires 
the lattice Mach number to be less than one. That is, in units where Ax = At = 1 , the 
speed of sound 

v s « y/dp /dp < 1. (17) 

In Lennard- Jones units v s ps \^50a/r for e* = 5 (see below) so we take At = O.lr so that 
v s « 0.7 Ax/ At. 

As discussed in Section IIV| the parameters determined from molecular simulations turn 
out to be qualitatively different from those commonly used in lattice Boltzmann simulations. 
As such, we had some problems with numerical stability using standard schemes. Stability 
was improved by using the predictor-corrector scheme |30J, rather than the standard Euler 
scheme. Stability can be further enhanced by iterating the corrector step a few times. This 
was found to be helpful in the initial steps, especially if a particularly poor initial state 
was used. In addition, the method for discretization of derivative operators, particularly 
Laplacian operators, made a significant difference. Including a mixture of derivatives along 
coordinate directions and those taken along the diagonal direction improved stability. 

To fully specify the model for the lattice Boltzmann algorithm, in addition to Eq. (|T2j) for 
the pressure tensor, we need explicit expressions for /i$ and po. These will be derived and 
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fit in later sections but for reference we list the complete expressions here: 

2A 

-i^V 2 $ - K p s,V 2 P) (18) 
Mo . If dA 2 \ 

Po = p W~ Ao + z{ p ^i- A2 ) x 



where 



[e0 + $ co ) 2 + e r ($ - $ co ) 2 - 2e r ei($ 2 + $ 



Z dp 



[ej($ + $ co ) - e r ($ - $ co ) - 2e r e z $] 



-P (k^p + i^V 2 $ + ^(V$) 2 ) , (19) 



e r = exp(-A 2 ($ - $ co ) 2 /(2pk B T)), 
e x = exp(-A 2 (<5> + <5> co ) 2 /{2pk B T)), 

Z = e r + t\ — e r ei. (20) 



III. MOLECULAR DYNAMICS MEASUREMENTS 

As seen in Fig. there are a number of parameters which affect the phase (mixed or 
separated) of the system. Much of the interesting behavior in binary fluids occurs when 
two phases coexist, so we will work with systems near the coexistence line. Coexisting 
phases have the same T and e*, but p and $ vary through the interface. Unless otherwise 
stated, simulations are performed at ksT/e = 1.1. The free energy functional Eq. 01 is then 
determined from simulations at e* = 5 and 6. 

The location of the coexistence line <3> co (p), a section of which is shown in Fig. El can be 
easily found by varying the average density and letting the system equilibrate (a process 
we accelerate using the Monte Carlo routines). To obtain further information about the 
free energy in the vicinity of an equilibrium state, we first study the linear response of the 
system to small perturbations. As we show in Section IIII Al this furnishes us with the 
second derivatives of the bulk free energy ip and the elastic constants in the equilibrium 
state. Section IIII Bl considers interfacial behavior. The surface tension is evaluated as a 
function of p and T, and the effect of capillary waves on the measured interface shape is 
examined. 
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FIG. 3: Portion of the coexistence curve for e* = 5(*) and e* = 6 (♦). The solid line is a quadratic 
fit to the data and the dashed line is from a fit to Flory-Huggins theory described in Appendix iBl 
Error bars are comparable to the symbol size. 

A. Linear Response of Equilibrium States 

We begin with an equilibrated single phase configuration from a molecular dynamics 
simulation. Thus the 'basis' state has no gradients of any kind. We then add a small 
force on atoms of type i, F% = SiCos(qx). The forces can be incorporated into the chemical 
potentials, Eq.(jSJ), as 

Si + 5 2 



F 



Hp 



2qk B T 

Si - s 2 



sm(qx) + 0(S 2 ), 
sm(qx) + 0(S 2 ), 



(21) 



where \i p and /x$ are given in Eq.J 

P = Po + Ps sin(qx) + 0(S 2 ) 
$ = $ + $ 5 sin(gx) + 0(S 2 



2qk B T 

If the Si are small, we expect a linear response. Then 



(22) 



where po and $o represent the undisturbed state. Plugging this into the previous equations 
and expanding about the equilibrium state gives 

+ K p q 2 jgfe + JW 
& + K*q 2 

where the derivatives are about the equilibrium state. By varying the S{ and wave vectors q 
we can determine the second derivatives of ip and the elastic constants K of the equilibrium 
state. 
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FIG. 4: Linear response coefficients as a function of wave vector squared. The slope of each linear fit 
(lines) gives an elastic constant and the intercept gives a second derivative of ip. The temperature 
is k B T/e = 1.1 and e* = 6. Average density is 0.83 (♦), 0.85 (*), 0.87 (■), or 0.9 m/a 3 (A). 
Statistical error bars are comparable to symbol sizes. 

To separate out the q dependence we use the following technique. For a given q we 
run two simulations, one with 5i — 82 — 0.2 and another with 8\ = —62 = 0.1. This 
gives us four equations to determine the four elements of the coefficient matrix (in reality, 
we have one extra equation due to the symmetry of the matrix). This process must be 
repeated for a number of different q in order to separate out the q dependence of each term. 
We obtain results for six g's simultaneously by apply perturbations at two different (and 
incommensurate) q's in each of the x, y, and z directions. 

The results for L can be seen in Fig. 0] The first thing to note is that the matrix 
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elements are indeed linear functions of q 2 . This indicates that the square gradient theory, 
assumed in Eq.flHJ), is capable of describing this system even at wavelengths as short as 2a. 
A surprising result is that both L pp and L p $ have a negative slope, implying that K p and 
Kp$, are both negative. This means that the square gradient theory becomes unstable to 
fluctuations at short wavelengths (~ la) where one of the eigenvalues of L crosses zero. As 
this distance is comparable to the molecular separation, fluctuations on shorter scales are 
unphysical. However, this sets a hard lower limit for the mesh spacing introduced in the 
lattice Boltzmann simulations. Using a mesh spacing shorter than the typical molecular 
separation to "improve accuracy" is not only pointless, but will be unstable jsil ]. 

It is also interesting to note that L pp is considerably larger (5 — 10 times) than the other 
components. This reflects the fact that it is harder to create fluctuations in density than 
fluctuations in concentrations. Conversely, a small change in density, such as the one seen in 
Fig-Eat an interface, may cost a non- negligible amount of free energy and thereby contribute 
significantly to the surface tension. We will discuss this point further in the next section. 

The second derivatives of the free energy and the elastic constants can be obtained from 
the linear fits in Fig. |3J Fig. El shows the second derivatives of if) as a function of p, and 
Fig. El shows the elastic constants plotted against p and/or <f>= $/p. Since values of ip and 
the ICs were only evaluated near the coexistence line, $ co (p), the dependence on density 
and concentration is intertwined. The detailed procedure for finding the functional form of 
ip is discussed in Sect. IIVI We focus here on parameterizing the elastic constants because 
they are important in understanding the interfacial tension results in the next section. Note 
that the variation of the K 7 s with p and $ was ignored in the chemical potentials of Eqs. (J7|) 
and (|H|). Appendix 1X1 explains why the variations discussed below can be ignored in these 
equations. 

To separate the $ and p dependence of the elastic constants we must consider points 
off the coexistence line. In principal this requires simulations on a two-dimensional grid. 
However, we find that K p appears to depend only on p, while K p $ and K$ depend only 
on the scaled variable = $/p. This is illustrated by including a point, p = 0.83m/<r 3 , 
$ = 0.79m/cr 3 , in Fig. El that is away from the coexistence line. This point collapses on 
the remaining data when K p is plotted against p and K p $, and are plotted against <fi. 
However, it lies far from the other data when the elastic constants are plotted against the 
other variable, as illustrated in the plot of K$> vs. p. Perhaps most surprising is that the 



14 



70 

CM 

E 

CO 

b 60 








(a) 

A. 
m 


50 










40 


— < — i—i — 1— 


, , 1 


i i — I— i- 




E^ 
b 

© 


★ 


X ★ 


★ 


(b) 


-15 


1 1 1 


i — i — i — I— 


-4 1 1 1 1 — 




10 








/ ' n 
/ * 


CM 
CO 

b 

"1 5 








' # / 

/ 7 

/ y 
y y 
* 

y 


CO 


^ ★ " " 


-* ' 


— S 
^ ♦ " 





0.84 0.86 0.88 0.9 0.92 



FIG. 5: The second derivatives of ip as a function of density along the coexistence line for (♦) 
e* = 6 and (*) e* = 5. Data points were obtained from the intercepts of the linear fits in Fig. 0] 
The second derivatives of the parameterization of ip for e* = 6 (solid line) and e* = 5 (dashed line) 
given in Section llVl and the Flory-Huggins free energy (dotted line) discussed in Appendix iBl for 
e* = 6 are also shown. 

results for different e* also collapse in Fig. |U1 (a), (b) and (d). Thus we can fit the data for 
the elastic constants with one set of parameters. 

Given the above results we take K p independent of $ and fit it to a linear function of 
p with the result given in Table By symmetry, K p <$> must be an odd function of $. The 
dashed line in Fig. Efb) corresponds to a fit of the data assuming K p § ~ whereas the solid 
line is a fit assuming K p ^> ~ 3 . The cubic function yields a better fit, giving the result 
listed in Table |l] This suggests that the common practice of neglecting terms involving K p <$, 
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FIG. 6: The elastic constants K p , K$>, and K p § as a function of the density p or order parameter 
(j) = $/ p. The points with the □ symbol are for e* = 5 and the ♦ and * symbols are for e* = 6. All 
systems lie on the coexistence line except the * point is off coexistence at p = 0.83, <3? = 0.79 m/a 3 . 
All elastic constants are in units of a 5 /m 2 . The lines represent fits given in Tableland described 
in the text. 

is justified in situations where is small, such as near the critical point, but is not justified 
in the more typical case considered here. 

Symmetry constrains Kq, to be an even function of <fi. The dotted line in Fig. Efd) comes 
from a fit of Kq, to a quadratic function of <ft. This clearly deviates systematically from the 
data which appears to be approaching a much flatter function of <ft as <fi — > 0. The solid and 
dashed lines, which both furnish good fits over the range of data, correspond respectively 
to the 8 and <ft 10 fits given in Table As we shall see in the next section, the value of 
Kq, has a significant impact on interface properties. However, these properties are most 
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(3.16 ±0.5-4) + (-6.08 ±0.57 4) P 




(-0.225 ± 0.011 ^)</» 3 




(0.286 ± 0.015 ^) + (0.517 ± 0.031 <p 8 




(0.319 ± 0.013 -4) + (0.524 ± 0.03 4) </> 10 



TABLE I: Elastic constants measured for e* = 5 and 6 and 0.82 m/a 3 < p < 0.925 m/er 3 . Care 
should be taken in extrapolating the elastic constants outside of the measured range. The two 
alternatives given for K§ fit equally well over the range of the linear response data (0.74 < <p < 
0.97). 

strongly dependent on the value of Kq> around = 0. As such, the fact that K<$> is nearly 
constant in this region is a highly desirable result. However, in order to have confidence in 
an extrapolation of Kq> from values of > 0.74 down to <p = requires that we examine the 
interfaces explicitly. 



B. Interface characterization 



We will examine the geometry shown in Figure |2] There are periodic boundary conditions 
in all directions and two flat (on average) interfaces normal to the x-axis. We have exam- 
ined different runs with average densities ranging from p = 0.82m/ a 3 to p = 0.925m/ a 3 , 
temperatures from ksT/e = 0.8 to ksT/e = 2.0, and e* from 5 to 8. Fig. [7| shows $ and p 
for a typical system at fc^T/e = 0.8. 

Most models of binary fluids in the literature assume incompressibility and therefore that 
p is constant. This is a reasonable assumption in the bulk but is violated at an interface, 
as illustrated in Fig. [7| The density dip in panel (b) is a ubiquitous feature of fluid-fluid 
interfaces noted recently in several publications j32,l3|, including some involving polymer 



fluids 



14. 



34j . The dip occurs because the energetically unfavorable 1-2 particle interactions 
are concentrated at the interface. The system can lower its free energy by decreasing the 
overall density at the interface as this reduces the number of 1-2 particle interactions. The 
size of the dip is determined by balancing this energy gain against the entropic cost from 
denying particles at the interface some volume. 

The surface tension associated with the interface is the most important characteristic for 
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FIG. 7: (a) Order parameter <3? and (b) density p along the x-direction of the simulation box, and 
averaged over the y and z directions. The temperature is T = 0.8 and e* = 5. (c) versus p for 
different interfaces with T = 1.1 and average density of 0.83 (•), 0.85 (♦), 0.87 (*), 0.9 (■), and 
0.925 (A). 

determining macroscopic behavior. The mechanical definition due to Kirkwood and Buff 
3o| relates the surface tension to the integral of the difference between the normal P± and 
parallel Py components of the pressure tensor across the interface. For a flat interface normal 
to the this can be written as 

7 = / [Pxxix) - (P yy (x) + P zz {x))/2] dx. (24) 

In the quiescent state, P xx is constant throughout the system (no flow implies d a P a p = 
and all quantities are constant in both the y and z directions). 
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Far from the interface the pressure is isotropic and one expects P xx = P yy = P zz . However, 
for small systems this expectation fails. For example for L y = L z = 8.2a and L x = 37.8a we 
found that P± — P\\ = 0.004 in a homogeneous system (i.e. without interfaces). Although this 
sounds small, when it is integrated over the whole system it yields (P± — P\\)L± = 0.15. This 
integrated stress difference is not related to the surface tension but still gives a significant 
contribution to Eq.(J2U). The effect appears to be due to the fact that the density- density 
correlation function has not decayed to zero at a separation of L y = L z but has by L x |l6[ . 
This appears to be a significant unrecognized error in molecular dynamics calculations of 
surface tensions. For the larger systems for which we present data, we find that (P± — P\\) < 
0.0002 and (P ± - P {] )L ± < 0.015. 

Fig. |S] shows measured surface tensions as a function of density, e*, and temperature. 
Normally, one expects the surface tension to decrease as you move towards the critical 
point. This is certainly the case as one decreases the density (Fig. Eta))- However, the 
surface tension increases as a function of temperature despite the fact that increasing the 
temperature moves the system towards the critical point. This effect has been noted before 
by a number of workers The surface tension peaks around ksT/e = 2 and then 

starts decreasing 3^. As one increases T, there is a stronger mixing at the interface that 
introduces more of the weaker 1 — 2 interactions and raises the potential energy. This leads 
to an increase in 7 until at higher temperatures the entropy contribution brings it back down 
again Q|. 

In addition to the surface tension, it is possible to define a characteristic width to the 
interface. The interface profile for $ from the molecular dynamics simulations (Fig. UZ[a)) 
can be fit to a tanh shape or equally well to an error function erf ((x — Xq)/ a/2£ 2 ). Similarly 
p can be fit to a constant minus a Gaussian, v4exp((:r — xo) 2 /(2£ 2 )). The resulting widths 
£ are shown in Fig. 

The stress profile through the interface also provides useful information. The surface 
tension is related to the stress through the kernel of Eq.()24|). 

r(x) = P xx {x) - {P yy {x) + P zz (x))/2, (25) 

We measure the local components of the stress tensor directly in the molecular dynamics 
simulations. A typical example is shown in Fig. IPOT a). The variation of T is also fit to a 
Gaussian to obtain a width. As shown in Fig. El this width is smaller than those obtained 
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FIG. 8: The surface tension 7 for e* = 5 (■*•) and 6 (♦) as a function of (a) average density at T = 1.1 
and (b) temperature at p = 0.85. Error bars represent systematic errors from {P xx — (Pyy + Pzz)/2) 
in the bulk regions. The dashed and solid lines in (a) are the surface tension computed using the 
lattice Boltzmann method for e* = 5 and 6, respectively, as discussed in Section The dotted 
line is the surface tension computed using the lattice Boltzmann method with the Flory-Huggins 
free energy discussed in Appendix El for e* = 6. 

from the interfacial profiles of the p and $. We will discuss this difference below. 

The interfacial stress profile is also related to gradients of the density and order parameter. 
Substituting the expression for the pressure tensor Eq.(JT2J) into Eq.(J25'|) gives, 

T(x) = k B T [K p {d x pf + K^d^f + 2K p *{d x p){d x §)} . (26) 

The numerical derivatives are computed by first doing a local quadratic fit to the data (p or 
$) over a range of five to seven points (data was collected in bins with width 0.29a) and then 
using the derivative of the local fitting function. Figure ITf¥ b) shows the contributions from 
the various terms in Eq. As is evident from the figure, the only significant contribution 
to the interfacial stress comes from the Kq>(d x §) 2 term. As a result, since the interface 
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FIG. 9: The interface width for e* = 6 as a function of average density calculated from fits of the 
order parameter <J> profile to an error function (♦) and the density profile to a constant minus a 
Gaussian (*). The width obtained from the stress profile T is also shown (■). Lines are only a 
guide to the eye. 



profile for <£> is well fit by the error function erf ((x — xq)/ a/2£ 2 ) we would expect the width 

obtained from a fit of the stress difference T(x) to a Gaussian to give a width l/y/2 « 0.71 
of that from the fit to $. The width we actually measure from T(x) is about 0.9 of that 
obtained from $ (see Fig. EJ). That is, the stress profile is wider than expected, a fact that is 
also obvious from Fig. 110( b). This may be due to a failure of the square gradient theory or 
may be a result of interface fluctuations not yet taken into account, namely capillary waves. 

The widths measured in molecular dynamics simulations are not the intrinsic values but 
widths broadened by thermal fluctuations. Capillary waves result in the measured time 
averaged width £ being larger than the intrinsic width £ by j3| 



where Ao is a short scale cutoff. Ao is expected to be proportional to the intrinsic interface 
width (fluctuations on length scales shorter than the interface width are not capillary waves). 
We take A = c£o where £o is the intrinsic width of the $ profile and c is a numerical constant. 
The bare width from the $ profile is used for £o since the width from T is a derived quantity 
in the square gradient theory. As all widths in the problem are proportional to the intrinsic 
width of the $ profile, changing the width used in defining A will change c but not A . 

Following Ref . jsfij] we verify Eq.(j2Zj) by plotting the measured width as a function of 
system size. This is shown in Fig. ITTTa) . Fitting the data to the expression £ 2 = ag + b§ In L, 





(27) 
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FIG. 10: (a) The interface stress difference T(x) = P xx (x) — (P yy (x) + P zz {x))/2 for a system with 
average density p = 0.83 and e* = 6. (b) The interface stress difference T (♦) and square gradient 
contributions K$(d x &) 2 (*), K p (d x p) 2 (□), and K p $,(d x p){d x &) (a). The elastic constants were 
taken from the linear response data (the <p 8 fit was used for K$). Data is for a system with average 
density 0.83 and e* = 6. Note that due to the periodic boundary conditions there are two interfaces 
in the system, giving two independent curves in (b). 

we find that b$ = ksT '/ '(2^7) to well within measurement errors, as expected from Eq. ([27)1 . 
The other parameter from the fit 

as = £0 -Mnc£ (28) 

could be used to obtain £0 if we knew c. Unfortunately that information is not available from 
this fit alone. In Ref. c = 13 was picked arbitrarily as the smallest number where Eq. (J28|) 
had a solution for all the systems considered, including monomers and polymers with up to 
30 monomers. There is, however, no reason to believe that c should be a universal number 
in this sense. 

If we assume that square gradient theory describes the interfaces, just as it works for 
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FIG. 11: (a) Dependence of the interfacial widths computed from the <J? profile (♦), p profile (*), 
and r (■) as a function of the system size parallel to the interface. The density is 0.85 and e* = 6 
for all systems. Lines show linear fits to Eq. 1271 (b) The intrinsic interface width as a function of 
average density for the order parameter $ profile from molecular dynamics (♦) and from lattice 
Boltzmann (line), (c) Parameter c relating the short scale cutoff Ao to £$ calculated using Eq.(J2HJ) 
(♦) and Eq. (|28|) (□). Data for e* = 6 is shown but similar results are obtained for e* = 5. 

linear response, then the value of c can be obtained from the difference in the widths of the 
T and $ profiles. In the square gradient theory the intrinsic widths should obey the relation 
£r = £i /2> where £r is the intrinsic width of Y and the intrinsic width of $. Based on 
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Eq.(|27|). the measured widths, and £r, of the $ and T profiles should be 



2tt7 V^*o 

Solving these equations simultaneously for c and £<j, gives the result shown in Fig lTTT b). 
Using Eq. (|28j) . we can also compute c from the finite size scaling if we use the intrinsic widths 
obtained from Eq. (j29j) . The agreement seen between these two methods in Fig. ITTT c) gives 
weight to the assumption that square gradient theory describes the interfaces accurately 

We can also assess the consistency of the values of K§ obtained from the linear response 
with estimates from interface data. To do this, we assume capillary waves broaden the 
interface in a Gaussian manner (i.e. the measured shape of the interface is a convolution of 
the intrinsic shape and a Gaussian distribution of width a/£ 2 — Q)- If we further assume 
that Kq> is constant and that the intrinsic line shape for $ is reasonably approximated by 
an error function (which is the case for the broadened shape directly measured), then 

^ = 0) J = k B TK* J* . (30) 

The quantity on the left is the directly measured (capillary broadened) stress difference 
divided by the gradient of $ at the point where $ crosses zero (center of the interface). The 
ratio of the length scales, £|/(\/2£r£$o)' is about 1.21 ± 0.05 for the interfaces measured. 
This gives a value of K$ of 0.29 ± 0.03<r 5 /m 2 . As T (and d x <&) is peaked at $ = and 
goes to zero as $ — > $ co this measure of K§ is dominated by the value of K<$> at $ = 0. In 
the previous section we determined that Kq> is essentially constant near $ = and that in 
this region it should have a value around 0.286 to 0.319 <r 5 /m 2 . Thus we can be reasonably 
happy that the elastic constants determined from linear response are consistent with those 
obtained from the interfaces. 



IV. FREE ENERGY PARAMETERIZATION 

In order to use the information from the molecular dynamics simulations in a mesoscopic 
model, we also need to parameterize the bulk free energy density ip. As we have direct in- 
formation about if} and its various derivatives from the interfacial stress and linear response 
data, this would not at first appear to be a difficult thing to do. However, if> is a function 
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of both p and $. The addition of the order parameter means that expressions for obtain- 
ing the free energy from the pressure found in standard references [l8| are not applicable. 
Fig. [7(c) shows that as one goes through an interface, $ can be reasonably approximated 
by a quadratic function of the local density p. As a result, distinguishing r/>'s dependence 
on $ from that on p using interfacial data alone is extremely difficult. Also, capillary waves 
change the shape of the interfacial profiles so that getting information about if) directly from 
the interface shape can be misleading. We therefore use the linear response data to do the 
fits and then compare the resulting surface properties as a test of the parameterization in 
Section El 

Before deciding on a particular functional form, it is worthwhile to examine what we 
would like to fit, and to prioritize what is most/least important. Our fit priorities are: 

1. Phase diagram/coexistence line ($ co (p)). 

2. Linear response of equilibrium phase (compressibility, etc). 

3. Surface tension. 

4. Interface widths. 

The rationale for this choice is fairly simple. Any model incapable of satisfying the first item 
has little, if any, useful properties. A model capable of reproducing the first two items will 
at least have a reasonable response to bulk forces. If a model fits the first three items then 
it should reproduce most macroscopic interfacial behavior. If in addition the fourth item 
can be fit, then we can reasonably expect it to describe a number of microscopic processes 
with macroscopic consequences, such as droplet coalescence or pinch-off. 

There are a number of free energy functionals commonly used to study fluid mixtures. 
We found Flory-Huggins theory has too few parameters to satisfy our first two requirements, 
while a Landau- Ginzburg expansion requires too many, resulting in ambiguous fits. In this 
section we discuss an alternative parameterization and relegate discussions of Flory-Huggins 
theory to Appendix [Bl 

To minimize degeneracies in the fitting procedure, we will try to parameterize the free 
energy in terms of quantities that are directly measured in the molecular dynamics simula- 
tions. The partition function we use is a sum of Gaussians centered at ±$ co plus a third 
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Gaussian centered at $ = whose amplitude is chosen so as to partially cancel the overlap 
of the other two Gaussians. This results in the free energy density: 

k B Ti) = A -pk B T\n[exp(-A 2 (<S>-<S> co ) 2 /(2pk B T)) 
+ exp(-A 2 ($ + $ co )7(2pfc B T)) 

-exp(-A 2 (<S> 2 + $ 2 co )/(pk B T))} (31) 

where Aq, A 2 , and $ co are functions of p to be determined. As we will show, Eq.(|31)l 
reproduces the linear response data to within statistical errors essentially by construction. 
For a single phase system on the coexistence line, but far from the critical point so <3> co >> 0, 
Eq.()31jl simplifies to a more familiar form 

1 



k B Ti) = A + - A 2 ($ 



$co) 2 . 



The parameters of the free energy functional can be determined directly from our molec- 
ular dynamics measurements of the equilibrium state. On the coexistence line ($ = ±$ co ) 
one can easily show that 



dp 2 

aim 



1 d 2 A + d 2 ip fd$ 



k B T dp 2 d$ 2 
d 2 ip (d® 

A 2 



dp 
1 — exp 



dp 



2Ao$l 



(32) 
(33) 
(34) 



k B T \ \ pk B T 

For states where $ co is significantly greater than zero the last expression reduces to 
d 2 i/j/d& 2 = A 2 . The equilibrium hydrostatic pressure is 

dA 



Po = P- 



dp 



-A . 



(35) 



By rearranging these equations it is also straightforward to show that along the coexistence 
line 

21 



pk B T 

d 2 ip 



dpo 
dp 

d$rn 



dp ' ■ 



d 2 i) 

dp 2 



( V 

\dpd$> J 

a* 2 



P- 



d 2 A 



o 



dp 2 



dpd<i> 



(36) 



(37) 



26 



0.84 



0.86 



0.88 



0.9 



0.92 



P 




0.84 



0.86 



0.88 



0.9 



P 



FIG. 12: Fits to (a) po and (b) dpo/dp which give the parameterization of Aq in Eq.(J3BJ))- In (a), 
molecular dynamics data is shown for e* = 5 (■) and e* = 6 (♦) and the solid and dashed lines 
correspond to the fits described in the main text. The dotted line in (a) (which overlaps with the 
solid line except at small p) corresponds to the fits used to determine V'o f° r the Flory-Huggins 
model described in Appendix [B] for e* = 6. In (b) only the e* = 6 data is shown. 

To use Eq. fl31|) in the lattice Boltzmann algorithm, we need to choose a functional form 
for the dependence of the parameters on p. As we have information about both values and 
derivatives (of $ co for instance), one could use a spline fit to match the measured values 
exactly. Alternatively, we can pick a functional form for the parameters and fit them over a 
range of p. For A , a convenient form is 



-4, 



o = ao + ciip In p + a 2 p 2 . 



(38) 



Then from Eq.fjHSJ), 



P0 = -O-O + CLlP + &2P 2 



(39) 
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and from Eq.(|37|) 



ldp d Z 1p \dpd$) 1 



ai - + 2a 2 . (40) 



p dp dp 2 ( d 2 ip\ p 

A simultaneous fit of po and dp^/dp, weighted by the statistical errors from the measurements 
of po an d the second derivatives used to obtain dpo/dp, is shown in FigIT2l The values of 
the fit parameters are given in Table ITT1 

Similarly, it is straightforward to do a fit of $ co to a quadratic function taking into account 
the information about d$ co /dp (Eq. (|37j) ): 

$co = b + b 1 (p-p ) + b 2 (p-po) 2 , (41) 

where po = 0.85 m/a 3 is a reference density chosen so as to make the statistical errors 
in the fit parameters given in Table |H] independent. This fit is shown in Fig. EJ The 
parameterization of the remaining quantity, A 2 , is done in two steps. First, Eq. (|34|) 
is numerically inverted to obtain A 2 at each point using the directly measured $ co and 
d 2 ip/d§ 2 . The resulting values are then fit to the function: 

A 2 = c exp (ci(p - po)) • (42) 

The fitted parameters are given in Table HP and the resulting fit to d 2 ip/d§ 2 is shown in 
Fig- Etc)- An exponential was chosen to fit A 2 rather than a quadratic or other polynomial for 
two reasons. First, it gave a better fit, and second, the quadratic did not remain monotonic 
over the full range of interest leading to unphysical effects when extrapolating outside the 
range of densities where linear response was measured. 

One could, in principle, add additional terms to the bulk free energy density that would 
have no impact on the coexistence curve or the linear response data fitted so far. Such a 
term would be zero and have zero first and second derivatives on the coexistence line. If, 
in addition, it was peaked at = it would affect the surface tension and interface width. 
As discussed in the next section, we do not need such terms here but they may be useful in 
other contexts. 



V. TESTS OF PARAMETERIZATION 



To test the above parameterization's description of interface properties we make use of 
lattice Boltzmann simulations. Lattice Boltzmann simulations were run for the same system 
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-9.5 ± 1.6 ^ 
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2.011 ± 0.060 -2^ 


2.827 ±0.007-^ 


Cl 


20.40 ±0.67^ 


21.70 ± 0.09 s£ 



TABLE II: Table of parameters for the free energy -0 in Eq. ()31|) . Parameters are defined in Eqs. 
(j55}) , (|4*T1) , and (|4*2*)) . As for the elastic constants given in Table QJ the data were measured for 
0.82 m/a 3 < p < 0.925 m/cr 3 and care should be taken in extrapolating outside of the measured 
range. 

size and compositions as the molecular dynamics simulations. We first verified that the 
lattice Boltzmann program correctly reproduced the bulk data. Surface tensions were then 
computed from the integral of the stress difference through the interface, Eq.()25|). just as we 
did for the molecular dynamics simulations. Figure |H1 shows the surface tension computed 
from the lattice Boltzmann program using the elastic constants and ip parameters from 
Table IHl Although the parameterization only used the linear response data, the surface 
tensions are very close to the molecular dynamics results, although they tend to be slightly 
higher. As discussed in Section IIII B[ the interfacial properties are dominated by the peak 
in the free energy and value of the elastic constants near = 0. All the fits to the data were 
for values of > 0.73, and it is remarkable that such an extrapolation works so well. 

Due to capillary waves we cannot directly compare the full interfacial profiles of the 
density and order parameter to the molecular dynamics. However, we can verify that the 
total mass deficit at the interface is correct. This can be calculated from the integral under 
the dip in plots like Fig. Efb) or from the increase in the density far from the interface that 
compensates for the dip. For example, consider a molecular dynamics simulation where the 
total average system density is p = 0.85m/ a 3 , k B T/e = 1.1, and with length L x = 74.36a. 
We find that the density far from the interface p = 0.8526 ± 0.0001m/a 3 is larger than the 
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FIG. 13: Order parameter profile on a cut going through the center of the drop. Note that the 
molecular dynamics profile has been smeared somewhat due to capillary waves and Brownian 
motion of the drop that are not present in the lattice Boltzmann simulation. 

system average to make up for the deficit at the interface. A lattice Boltzmann simulation 
with the same system size and total average density shows exactly the same mass increase 
in the bulk. 

Figure UTT a) shows that there is also reasonable agreement between the intrinsic widths 
computed via the lattice Boltzmann and the molecular dynamics. Although in this case the 
intrinsic widths measured from the lattice Boltzmann tend to be a bit smaller. 

To test our model parameterization in a different geometry we examined a cylindrical 
drop. In three dimensions, a cylinder of fluid is unstable to spherical droplet formation. We 
avoided this instability by making the radius R of the droplet bigger than the system size 
along the axis of the cylinder (y-axis). Specifically, R = 24.5a, L y = 16a, and L x and L z 
were about 100a. For a density of p = 0.85m/ a 3 this required 147456 molecules and 10 7 time 
steps were needed to get good statistics, which is reasonably large for a molecular dynamics 
simulation. The only significant difference between the molecular dynamics simulation and 
the lattice Boltzmann method is the presence of thermal noise. In the molecular dynamics 
simulation this causes the drop to undergo Brownian motion. In order to do a meaningful 
comparison we limited this effect by periodically shifting the system so that the center of 
mass of the drop remained at x = z = 0. 

Figure EH compares profiles of the order parameter through the center of the drop from 
molecular dynamics and lattice Boltzmann simulations. The molecular dynamics profile is 
broadened slightly by capillary waves and Brownian motion. Neither effect is present in the 
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FIG. 14: Laplace pressure between the inside and outside of a cylindrical fluid "drop", versus 
the curvature 1/R, where R is the drop radius. Results from lattice Boltzmann simulations (*) 
and molecular dynamics simulations (□) are consistent with each other and with the straight line 
prediction from continuum theory Eg 1431 

lattice Boltzmann simulation. The two methods both show that <3> is lower by 0.01m/ a 3 than 
the equilibrium bulk value in the center of the drop. The reason for the drop in $ is that 
the interfacial curvature produces a pressure difference Ap between the inside and outside 
called the Laplace pressure. This also increases the density in the drop. For R = 24.5a the 
density difference from lattice Boltzmann simulations is 0.00867m/cr 3 , which agrees with 
the value of 0.009 ± 0.002m/cr 3 from molecular dynamics. For a macroscopic cylinder, the 
pressure difference between the inside and outside of the drop should scale with the radius 
as [37 1 

Ap = j/R, (43) 

where 7 is the surface tension measured in the molecular dynamics simulation of a flat 
interface. As can be seen in Figure El both the molecular dynamics simulation and the 
lattice Boltzmann results follow this relationship very well. 



VI. MODEL COMPARISONS 



It is worthwhile to compare our findings to commonly used assumptions in mesoscale 
modeling of binary fluids. We first consider the assumption of incompressibility that is often 
used for both simple and binary fluids. For bulk fluids, including those examined here, this is 
a very reasonable assumption as d 2 i/j/dp 2 is very large compared to other second derivatives 
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of the free energy (cf. Section UlI A|) . However, due to the fact that K p is negative, the fluid 
can be very soft on the short length scales characteristic of interface widths. This invalidates 
the assumption of incompressibility in the interfacial region. 

As noted in Section 111 Bl the density drops in the interfacial region to reduce unfavorable 
1 — 2 molecular interactions and hence the free energy. The quantitative impact of density 
changes can be seen by considering the change in the free energy barrier (Eq. I3~T|) in the 
interfacial region. Due to the change in density in the interface, we find that A2 drops down 
to around 50 percent of its value in the bulk. This reduction in A 2 spreads the Gaussian 
functions in Eq. (j31j) thereby reducing the barrier between the ±$ co states. There is a 
corresponding reduction in the total interfacial tension which is given by the integral of T(x) 
through the interface (Eqs. f2~il and 12*5]) . Using the Euler-Lagrange equations (Eqs. |7|and|SJ) 
one can show that T(x) is directly related to the free energy barrier: T(x) = ksTip — A$. To 
see the impact of density changes in the interface we can compare the integral of T(x) using 
the $ profile from the lattice Boltzmann simulation combined with either a p-dependent 
value of A2 or a constant bulk value of A%. We find that including density variations reduces 
the surface tension by a factor of 2 — 4 relative to incompressible models. Thus the density 
drop at the interface is not a negligible effect for any quantitative study involving interfaces. 

As mentioned before, the Flory-Huggins theory has too few parameters to obtain a precise 
fit to both the coexistence curve and the linear response. As discussed in Appendix |Bj we 
can obtain the Flory-Huggins \ parameter from a fit to the coexisence curve shown in 
Fig. El The resulting linear response seen in Fig. while not in quantitative agreement 
with MD results, should be adequate for qualitative work. The Flory-Huggins prediction 
for the surface tensions seen in Fig. |H] is actually remarkably good. It is important to 
note however, that we have allowed \ to vary with the local density. As a result, \ drops 
significantly in the interfacial region due to the interfacial density drop. This reduces the 
surface tension by a large factor compared to incompressible models, just as the reduction in 
A2 discussed above did. As most implementations of Flory-Huggins models for simulations 
assume incompressibility, such models will have surface tensions that are a factor of 2 — 4 
too large. 

It is commonly found in lattice Boltzmann simulations that a stationary droplet i 
cent conditions will develop a flow field around it similar to that shown in Fig. fTST a) 
It has recently been pointed out that these spurious velocities arise due to discretization er- 
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FIG. 15: $ field (shading) and velocity field (vectors) for one quadrant of a stationary cylindrical 
fluid "drop". In (a), a standard model from the literature [3| was used and in (b), the model we 
have matched to the molecular dynamics data. In units where dx = dt = 1.0, the largest velocity 
vector is 1.2 x 10 -4 in (a) and 1.5 x 10~ 6 in (b). Note that the velocities in (b) would not be visible 
if we had used the same scale as that used in (a). The radius of the drop in (b) is slightly smaller 
but this should only increase spurious velocities compared to (a). 
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rors 38] which drive the spurious currents. For our model, fitted to the molecular dynamics 
data, we can estimate the truncation error from the discretization to be ~ 10~ 6 and the ve- 
locities seen at the interface in Fig. lToT b) are indeed of this magnitude (in lattice units). The 
discretization error for the standard model from the literature used to produce Fig. IToT a) 
should also be ~ 10~ 6 however the spurious velocities observed are nearly 100 times greater. 
Further, the flow field for the standard model includes considerable vorticity and significant 
flows far from the interface itself. This suggests that discretization errors are driving an 
unstable (and possibly unphysical) mode of the standard system leading to much larger 
spurious velocities. In our model, this mode is absent and discretization errors are not 
exaggerated. 

One possible origin for unphysical modes in Fig. IToT a) and similar models is related to 
the time scales used. In situations where the model parameters have not been measured so 
that the length and time scales are uncertain, the viscosity and diffusion constant may be 
set for convenience. Normally this would not affect equilibrium properties, but if unphysical 
choices are made, unstable modes like the ones being driven by the spurious velocities may 
be set up. For instance, unstable modes may appear if the diffusion constant is set so large 
that material can diffuse faster than it can flow. To ensure that this is not the case in our 
model, we have used a viscosity and diffusion constant that were measured in molecular 
dynamics simulations for similar fluids, 7] « 3er/cx j3| and D = 0.1 a 2 /r Q|. 



VII. SUMMARY AND CONCLUSIONS 



This paper presents detailed molecular dynamics simulation results for a binary mixture 
of simple fluids and uses them to construct a square gradient theory that can be used for 
realistic mesoscale modeling. The MD simulations examine the linear response, interfacial 
tension, and interfacial width as a function of density, temperature, and the repulsion be- 
tween different species. Two remarkable conclusions arise from the linear response results 
(Sec. IIII A|) . The first is that a square gradient theory is capable of quantitatively describing 
the response down to wavelengths that are comparable to the molecular spacing (< 2a). This 
implies that mesoscale models may have a much wider range of applicability than might be 
expected. The second surprising result is that the elastic constant for density changes, K p , is 
negative. As a result, the system is more susceptible to fluctuations at shorter wavelengths. 
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Indeed, it is only stabilized at small scales by atomic discreteness [31] . This prevents density 
fluctuations on length scales less than a where the total response coefficient L pp becomes 
negative. 

Studies of interfacial properties (Sec. 1111 Bj) show that the common assumption of incom- 
pressibility is not valid. This is related to the observation that K p < 0, which lowers the 
free energy cost of localized density fluctuations. Although the density change is small, it 
can reduce the interfacial tension by a factor of 2-4. 

The density, order parameter and surface stress were evaluated as a function of position 
normal to the interface, and used to determine interfacial widths. The variation in width 
with system size is consistent with broadening by thermal capillary waves. Comparing the 
scaling of widths from the stress and order parameter, allows all the parameters of the 
capillary model to be determined independently. 

Fits to linear response about states on the coexistence curve showed an impressive ability 
to predict interfacial properties. Predicted values of the surface tension (Fig. (HJ) and the 
density deficit at the interface (Sec. EJ are nearly within the statistical error bars of the MD 
results. This is particularly surprising given the large change in order parameter through 
the interface and small interface width (Fig. Ej). The main discrepancy between the MD 
results and mesoscale simulations is that the latter do not include interface broadening by 
capillary waves. 

Many LB models have been found to produce spurious velocities around curved interfaces. 
While some discretization error is expected, it appears that these errors are amplified when 
the time scales in the LB model are chosen arbitrarily or for computational convenience. 
Using time scales derived from MD simulations prevents unphysical choices that, for example, 
allow material to diffuse more rapidly than it flows. Fig. IT51 shows that using MD parameters 
can reduce spurious velocities by around three orders of magnitude. 

The final square gradient theory (Table has a simple analytic form and provides 
an excellent fit to all MD results for phase coexistence, linear response, and interfacial 
properties over a wide range of densities. This is particularly important for future studies of 
nonequilibrium phenomena such as pinchoff or contact line motion. Dynamic processes will 
lead to variations in local density and order parameter that will in turn lead to variations in 
local interfacial stress. These variations will have important implications for the dynamics, 
producing Marangoni-like effects or changing the wavelength dependent response of the 
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interface. Our results for the influence of density changes on interfacial tension indicate that 
these effects may be quite large. While simpler theories (e.g. Appendix [Bj) may be able 
to describe equilibrium configurations, they do not include these important dependences 
on local density. More complex theories that are not guided by MD results are unlikely to 
include important effects such as a negative K p . It will be interesting to explore the dynamic 
consequences of such effects in future work. 
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APPENDIX A: ELASTIC CONSTANT VARIATIONS 

Frequently the free energy is parameterized as a function of the individual species densities 
Pi and pi- The free energy can then be expressed as 

T = k B T J dr |^ + I^ n (Vpi) 2 + 

\k 22 (V P2 ) 2 + K 12 V Pl ■ Vp 2 j . (AI) 

The resulting elastic constants are linearly related to those in Eq.(J3J): Ku = K p +Kg,+2K p g>, 
K n = K p + Kg,- 2K p g>, and K 12 = K p - Kg,. 

Another order parameter that is often used, the relative concentration, is = $/p. The 
free energy becomes 

T = k B T J dr jv + \k p {Vpf+ 

^ (V0) 2 + k p<t> Vp ■ J . (A2) 

These elastic constants have a more complex mapping to those in the main text, 

K p = k p + k^ 2 /p A - 2k p ^/p 2 , (A3) 
K* = k^/p 2 , (A4) 
K p g> = k p4> /p- k^/p 3 . (A5) 
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It is obvious from these relations that the ICs and fc's can not both be independent of p 
or <3>. There can be advantages to using either or $ in different physical situations. It 
turns out that for linear response measurements and for the lattice Boltzmann algorithm it 
is convenient to work with $. 

As the elastic constants can vary as a function of $ and p this should, in principle, 
be taken into account in Eq.fjHJ). For the elastic constants used in the main text, the full 
Euler-Lagrange equations are 



ldK r 



f p -K p ^p-K p ^--- 



(Vp) 2 + 



1 dK«, dK, 
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2 dp 

dip 



3$ 



(V$) 2 -^Vp- v$, 



ldK* 



(V$) 5 



1 dK n dK, 



2 <9$ 



dp 



(Vp) 5 



dp 



-V$ • Vp. 



(A6) 



(A7) 



Clearly it would be desirable if some of these terms, especially those involving variations of 
the elastic constants, were negligible. For the linear response regime it is straightforward to 
show that these additional terms are order S 2 (cf. Equation (J2HJ))- However, we must still 
investigate their importance for interfaces. 

Figure shows various contributions from gradients to the potentials of Eq. (|A6|) and 
(jA7j) in a typical interface. As can be seen, the terms that dominate have been kept in 
Eqs. and JHJ). The remaining terms are quite small, or essentially zero, thus justifying 
ignoring them in the main text. Note that some terms that are quite small in interfaces, 
such as K p ^V 2 p in p$, are required when looking at the linear response where Vp and V$ 
are of comparable magnitude. 



APPENDIX B: FITS TO FLORY-HUGGINS FREE ENERGY 

There are a number of free energy functionals in common use to study fluid mixtures. 
For polymer mixtures, the Flory-Huggins free energy density is most commonly used. For 
monomers the free energy density divided by k^T is 

^ = Piln(^)+p 2 ln(^ +X ^ (Bl) 
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FIG. 16: The individual gradient term contributions to the potentials (a) /x$ (Ea. (|A7|) ) and 
(b) p p (Ea. (|A6|) ) through two interfaces in a system at average density p = 0.85 with e* = 6. 
For 11$, the term if$V 2 3> (♦) dominates, T;(dK$/d<fr)(X7&) 2 (■) is small but visibly different 
from zero, and \(dK p /d<S>)(V p) 2 (*), (dK^/dpXVp) 2 (□), and (dK$/dp)V$ • Vp (a) are effec- 
tively zero, (b) For fi p the significant terms are K p V 2 p (♦), K p $V 2 ® (*) and (9/sT p $/(9$)(V$) 2 
(□). The ^(9X$/9 / o)(V < I ) ) 2 (a) term is small but visibly different from zero. The other terms, 
±(dK p /dp)(Vp) 2 (■) and (dK p /d<Z>)Vp • V$ (A) are effectively zero. 

where x is the only free parameter. This is just a modified entropy of mixing and, in order 
to obtain the bulk pressure, one must add to this an additional function ip of p alone, 



(B2) 



If one were to follow the spirit of the derivation of the Flory-Huggins model, i/jq should be 
determined primarily from the entropy of an ideal gas plus some quadratic terms to correct 
for energy interactions. In practice it is unrealistic to expect such a construction to work. 
We shall use the same form for ip as we used for A in Eq. (J38j) . There is also some ambiguity 
in the definition of the \ term in the Flory-Huggins free energy. Some authors use a slightly 
different term, XPip2/{vp 2 ) where v is a reference volume 40]. If we allow % to be a function 
of p then both terms are equivalent but the meaning of \ will be slightly different. 
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We obtain \ by fitting $ co and d$ co /dp as a function of p. In the bulk states p$ is zero 
and there are no gradients so that Eq.fjHJ) requires that on the coexistence line 

l + $/p 



= ®± = djjFH = l ln 
<9$ <9$ 2 



1 $ 

(B3) 

2 p 



.l-<&/p_ 

Eq. (}3Tj) also holds, and if one evaluates these derivatives for the Flory-Huggins free energy 
one can obtain the relation 

1 1 \ d$ r 



l + $/p l-$/p/ \p dp 

where d$ co / dp is evaluated using Eq. (J3T|) and the measured values of the second derivatives. 
If we take x to be a quadratic function of density, 

X = Xo + XiP + X2P 2 , (B5) 

then we can do a simultaneous fit to these two equations using a straightforward weighted 
linear regression (there should not be any conflict between them as we are just fitting $ co (p) 
and the derivative of this function d$ co /dp). The weights are computed from the statistical 
errors of the measurements of $ co and the second derivatives. We use standard methods to 
find the statistical errors of the derived quantities. The resulting fits are shown in Fig. El 
and the parameters are given in Table IIHI 

Next we fit the bulk pressure po to obtain ip . Unlike the fits of the equilibrium p in the 
main text, which involved only A , there is now a contribution from the $ dependent terms 
in the free energy. In particular, 

Po = P-q^ - ^0 + ^k B T{p - $ )— . (B6) 

One could, in principle, also use the information about dpo/dp that we used in fitting Aq for 
the single phase free energy. However, the additional terms make a combined fit extremely 
messy and of dubious value. One still needs a parameterization of i/jq i n terms of p and we 
use the same functional form as was used for A (although the numerical values for aoo etc. 
will, or course be different). The resulting fit to po is shown as a dashed line in Fig. [T21 and 
the parameters are given in Table 11111 

As all parameters in the Flory-Huggins free energy are now determined, we can now com- 
pare the quantities not explicitly used in the fits. Fig. |3] shows the second derivatives of the 
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FIG. 17: x parameter (a) obtained from the fits to Q co {p) (Fig. |3J) and d$ co /dp (b). Solid and 
dotted lines correspond to fits for e* = 6 and 5 respectively. Molecular dynamics data in (b) are 
for e* = 5 (*) and 6 (♦). 

bulk free energy as measured from linear response in the molecular dynamics simulations 
and from the fits. The Flory-Huggins theory overestimates the concavity of the free energy 
minima. Fig. |H] shows the surface tension derived from a lattice Boltzmann implementa- 
tion of the fit to the Flory-Huggins theory. The agreement is remarkably good. However, 
Flory-Huggins normally assumes incompressiblity. If we had assumed that the density, and 
therefore x was constant, the surface tension would be much too large. 

It is also worthwhile to compare the values of % obtained here to other methods of 
estimating \. For long polymers at p = 0.85m/ a 3 , Grest and coworkers |l3( have estimated 
X in two ways. Using a so-called one-fluid approximation, they estimate x ~ 0.76e*e/fceT. 
Using an incompressible random phase approximation to evaluate the static structure factor, 
they obtain a larger value of x ~ 1.0e*e//ceT. Our results correspond to a somewhat smaller 
value of about 0.54e*e//csT, which is not surprising given that we consider simple monomers. 
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e* 


5 


6 


«o 


-31.51 ± 2.9 


-26.50 ±6.3-^ 

(TT 


a\ 


-92.33 ±6.7 4 


-80.18 ±14 4 


a2 


77.0 ± 3.8 -4r 


69.7 ± 8.3 -4r 


Xo 


10.15 ±4.7 


10.68 ±5.0 


Xi 


-28.47 ± 10.4^ 

m 


-29.78 ± 11.21 ^ 

m 


X2 


23.42 ±5.7 -4 


24.52 ± 6.32 -4 



TABLE III: Table of parameters for the Flory-Huggins free energy given in Eq.s (|B1|) and (|B2f) . 
Parameters are denned in Eqs. (|38|) and ()B5|) . As for the parameters given in the main text, the 
data were measured for 0.82 m/a 3 < p < 0.925 m/cr 3 and care should be taken in extrapolating 
outside of the measured range. 
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